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CHAPTER 1 
Introduction - Qualitative Review 
1.1 Motivation 

Interest in the strength of naturally occuring ice and 
the ultimate size to which marine structures must be designed 
has become a major topic as the last boundaries of easily 
accessible petroleum resources are used. The extension of 
engineering practice into the world's colder regions is 
subsequently an effort that is intimately related to this 
perceived need. 

To satisfy a curiosity regarding this area of offshore 
design, the author set forth to discover a rational approach 
to establishing a meaningful standard. It was hoped that 
an appropriate theory, perhaps based on statistical patterns, 
was available that was similar to the many comprehensive 
theories regarding wind, wave, and seismic loading. Unfor- 
tunately, no comprehensive or universally accepted theory 
was found. However, many prominant investigators, including 
Assur [1 ][2 ][ 3], Bercha [4][5][6], Hirayama et a1. [14] 
[15], and Korzhavin [19][20], have proposed numerous models 
and approaches to the calculation of ice strength and the 
quantitative prediction of necessary design size. Some of 
these models will be reviewed in Subsequent chapters. 

Consequently, the motivation for this thesis is to 


recognize and segregate an accurate and useful theory for 
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the prediction of ice pressures, 
1.2 Scope 
The purpose of this investigation is to analyze the 
ice-structural interaction models presently proposed for tne 
case of a vertical circular pile. The circular pile was 
chosen primarily because it is the most useful geometric 
form in areas of high tidal variance. Other authors, 
notably Tryde [30] and Bercha [5], have examined the other 
promising forms, the inclined wedge and the cone. | 
The investigation will take the form of the following 
steps: | 
1) Qualitatively compare the available analytical 
models. 
2) Either improve upon the most promising model or 
construct a new model. 
3) Analyze sensitivity to some important parameters. 
4) Analyze some important ice properties affecting 
the circular pile interaction problem. 
1.3 Structural Failure 
The earliest designs contemplating an ice environment 
were highway and railroad bridges in North America and the 
Soviet Union. Some authors, notably Korzhavin [19] and 


Watts [32], have summarized structural collapses known to 
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be due to ice forces. These failures have a common thread; 
they were all caused by unusual ice jams which, coupled 

with an early thaw, produced unplanned loadings. Insights 
gained from a careful reading of these descriptions indicate: 

- very few failures actually occurred. This evidence 

implies that early engineers had either very accurate 
analytical skills (hardly likely) or that early 
designs were overly conservative. 

=- some failures were by gradual ice-erosion, 

= when failure did occur, it was due to an unusual 

meteorological complication or no allowance for 
ice at all. 

No major marine structural failures, apart from 
ice-bound ships, were found in the literature. Perhaps this is 
because this hostile environment has yet to be seriously 
challenged. Ice damage on coastline facilities, breakwaters, 
docking facilities, and lights, have been recorded. [26] 

1.4 Application 

Applications of ice engineering extend not only to the 
offshore industry but'over a wide range of related problems. 
The knowledge gained through the solution of marine ice 
uncertainty can apply to the following range of actual 
proposals: 

= over ice transportation 

- ice airfields 


=- submarine surfacing through ice 
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= jcebreaker technology - navigation 

- ice island construction for use as a moblie platform 

= iceberg towing 

The offshore oil industry has perhaps the greatest 
motivation in subduing arctic conditions, Although 
experience to date is not extensive, the few structures 
presently deployed emphasize the continuing need for more 
accurate analysis. 

Present North American interest is concentrated in 
Cook Inlet, Alaska, [26], the North Slope, and the Beaufort 
Sea [9]. Arctic conditions are most severe, es it will 
be shown in the next few sections, and there are designs 
under construction which do not even consider a passive 
survival of a structure, These proposals include a vibra- 
tory icebreaker motion [7], ice melting by internal heating 
[17], air cushion mobile platform [22], construction of 
protective soil berms [9], or even mobile platforms that 
only drill in optimum conditions, 

Scandinavian, Soviet, and Argentinian oil interests 
are not represented in the literature. 
1,5 Historical Treatment 

Analytical treatment of ice pressures and winter loading 
was first presented in fresh water civil works, especially 
dams and bridge piers. Early investigators faced the problems 
by correctly deducing that no matter what the loading condi- 


tion or geometric possibilities, ice strength could not exceed 
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its crushing strength. Thus a limiting value of 6% was 


imposed into the following relationship: 


Fax = bh Oca (1.1) 


b represents pier width and Gu. crush strength of.ice 


አ is force exerted on structure 


A special study of the A.S.C.E. in 1931 (Committee on 
Power Division) recommended a crush strength of 400 psi [32]. 
A review of Russian, Polish, Canadian, and U.S. military 
design codes was conducted by Watts [32], and it reveals 
that most design criteria are simply a refined empirical 
version of (1.1). 

Design criteria for offshore structures in the much 
more quantitatively unknown sea environment do not simplis- 
tically prescribe this relation. The paragraph relating to 
design for sea ice forces for fixed offshore structures from 
Det Norske Veritas 1974 Rules states [10]: 

B401. For structures intended to be installed in areas 

where ice hazards may exist, relevant statistical data 

for the area in question are to be submitted. The ice 
conditions are to be described with particular attention 

DOS 

- concentration and distribution of ice 

- types of ice (ice floes, ice ridges, rafted ice, etc.) 


- mechanical properties of ice 
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~ mean thickness of ice floes 

- average drifting speed and direction of ice 

- probability of encountering icebergs 

- tidal range 

This rather nebulous requirement is in stark contrast 
to the curves, tables, formulae, and data presented for other 
force loadings. It is here that the lack of analysis in 
the state of the art is most keenly felt. 

To increase the body of knowledge presently available 
in the science regarding ice pressure loadings, numerous 
investigators in the last decade have proposed many alternate 
procedures, The more promising approaches will be analyzed 
in & subsequent chapter. 

The proposed solutions to the ice pressure problem 
fall basically into four broad categories: 

l) Classical elastic thin plate theory 

2) Finite element or finite difference numerical solutions 

3) Empirical formulation 

4) Model basin testing 

Each approach has merits as well as severe, restrictive 
limitations, For instance, the classical theory assumes 
thin plate elastic behavior (plane stress or plane strain, 
isotropic, homogeneous, elastic) which is most certainly not 
the complete behavior of sea ice, Yet the solution can yield 


a fair degree of accuracy within certain limiting conditions 
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and provide a closed form expression which is easy to use 
and easily reproducible, 
1.6 Importance of Ice Forces 

To insure that the reader is qualitatively aware of the 
enormous magnitude that ice pressures can exert, a simple 
model is presented in Figure 1-1, 

A monopod platform with a radius at ice of 8 feet is 
subjected to an ice load driven by wind or current shear, If 
the extent of the floe is great enough, it is obvious that 
forces well in excess of local ice failure can be generated, 
Using a Oca of 200 psi and formula (1.1), adjusted for 


projected frontal area, the ice force can be: 
Lan E ES bh = OIO Tons 


This approximately 700 tons is concentrated on a three 
foot slice of the leg, presenting severe structural consid- 
erations. If this structure were located in 80 feet of 
water, the overturning moment could be 20,700 foot-tons, a 
serious concern to foundation and soil resistance as well as 
flexural strength of the leg. 

Perhaps an even more serious consideration which 
has recently been addressed [21] [25] [26] is the fact 
that this force can thrust at a frequency close to one hertz, 
Apart from resonant considerations, the fatigue strengths 
of the cold metal and cold weldments become of serious 


concerne 





XXX 
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Ice interaction with monopod type structure 
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Peyton [26] has summarized some important requirements 
for the ice force design of marine structures. These include: 

- Horizontal forces due to crushing, bending, buckling, 

or shearing of the ice sheet (quasi-static) 

- Impact forces (dynamic) 

= Weight of attached sea ice at low tide 

e Buoyant lift of above at high tide 

- Vertical components of sheet failure 

- Diaphragm bending during water level change 

- Ice accretion by spray 

- Thermal expansion of ice in joints 

=- Rubble in structural framing carried by ice 

e Abrasion by moving ice 

It can be seen that in a cold environment, ice forces 
can be large and obvious as well as an insidious problem. 
In the next section the enormity of the design equation 
and the limits of present knowledge are presented. 
1.7 Ice and the Enormity of the Design Problem 

Ice is a naturally occurring material which forms in a 
convenient plate rege located at the surface of the liquid- 
air interface, Ice formed on a small fresh water pond can 
be clear, homogeneous, constant in thickness and have a sen- 
Sible temperature distribution. In this case the engineer 
could possibly apply some of the more well-known classical 
plate solutions to find the ice behavior and predict with a 
fair degree of confidence any desired solution to his prob- 


len, 
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Sea ice is not that well behaved, however. Consider 


the situation shown in Figure l-2 which depicts a profile 


of arctic coastal ice, A bottom founded structure could 


be located anywhere in the picture. 


Figure 1-3 presents a careful summary of the unknown 


parameters affecting the design equation. A high degree 


of interdependence is characteristic, 


Figure 1-4 presents four logical paths to the final 


design. 


All have been used in successful designs. Again, 


each path has its advantages and disadvantages. 


1) 


2) 


3) 


From known environmental conditions - this is 
perhaps the ideal situation; a successful design 
could possibly be created anywhere at minimal 

cost, once perfected. Of course, the lack of 
knowledge regarding the environment at each locale, 
the rudimentary state of analytical tools, and the 
lack of correlation with known quantities makes 

this procedure years away, and perhaps unattainable. 
In situ - This will provide better data, but still 
suffers from analytical model imperfection, as above. 
In addition, the gathering of this data can be 
costly. 

Near full scale =- a very practical approach, but 

it too can be costly, nearly as much as the final 
product in some cases [25]. It does offer the 


advantage of risking a small, inconsequential failure 
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Typical winter features, Beaufort Sea 
(after Croasdale [9] ) 
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over a larger, more serious accident. 

4) Previous experience - This approach would yield 
the highest confidence level, but, unfortunately, 
little has been done and some of what has is 
proprietary information. Here again an overly 
conservative design would propagate other overly 
conservative designs. 

Without presenting a complete review of the parameters 
in Figure 1-3 or the literature on Figure 1-4, both of which 
are beyond the scope of this thesis, it is hoped that an 
appreciation of the enormity of the design effort involving 


ice can be gained. 
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CHAPTER 2 


Summary of Some Sea Ice Material Properties 


2.1 Purpose of this Chapter 
The purpose of this chapter is to present a brief. 


outline of some important sea ice parameters and properties 
for use in developing the analytical model presented later 
in this thesis. For a thorough treatment of ice properties, 
the interested reader is referred to references [12], [13], 
[18], and [33]. 
2.2 Problem Formulation 

Sea ice can be considered to be sensitive to temperature, 
salinity, crystal orientation, and seasonal variation. It 
can be considered to behave elastically for a short load 
duration and a moderate load magnitude. Katona and Vaudrey 
[13] describe an approach to an appropriate field theory 
which is sufficiently general to apply to the case of a 
circular pile - ice sheet interaction problem. The basic 
parameters are presented in Figure 2-1, 

Katona and Vaudrey further reduce the list of material 
parameters (category I of Figure 2-1) to a functional form 


which contains only two parameters: 


т 
/ / 
Qu 7 Ria Lea (t) x] + a (2.1) 
where Oi; is the stress tensor at spatial point x 


and time + 
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Parameters to define an appropriate field theory 
(after Katona and Vaudrey nal ) 
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m is à symmetric fourth order tensor functional 

allowing complete anisotropic formulation 

E. (f)is the strain tensor history O<t<t 

x is the spatial point included for any variation 

of other parameters 

Ë current time 
By exploiting knowledge of sea ice behavior, (2.1) can be 
reduced to a viscoelastic representation by restricting 
load magnitude and load duration. Further restriction 
results in an elastic formulation, The resultant material 
behavior can be characterized as falling somewhere into 
Figure 2-2, The quantitative limits of the regions have 
not yet been experimentally determined. 

2.3 Failure 

Category III of figure 2-1 requires knowledge of the 
mode of failure of ice surrounding a circular pile. Most 
investigators work with an assumption of compression failure 
of the ice immediately preceding the pile. However, since 
this failure is quite large with respect to shear or tensile 
failure in simple tests, it is conceivable that tension or 
Shear failure may precede compressive failure, 

Additionally, if it is assumed that the greatest level 
of force on the pile will be produced at the point of rupture, 
cracking and large deformation behavior need only be applied 
for postfailure behavior, 

Failure by limiting strain is possible also at levels of 


force below elastic limits and at larger time durations. 








load 
magnitude 






nonlinear viscoplastic 


constant temperature boundary 


effect of increased temverature 


vıscoelastic 


load duration 


FIGURE 2-2 


Effect on domain boundaries by load parameters at a 


constant temperature (after Katona and Vaudrey E ) 
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2.4 Material Properties 
There is a large gap in the body of knowledge regarding 


the material properties of ice and nowhere is there a con- 
sistent and universally recognized presentation of values. 
Therefore, for purposes of construction of the model and 
presentation of results of this thesis, non-dimensional 
quantities will be used wherever possible. Some representative 


quantities which may be applied for clarity: 


Quantity Typical Value 

Young's modulus E per 300,000 psi 
Poisson's ratio P «33333 
Shear modulus zi 112,500 psi 
Viscous constant X ger 1.2788 x 10? (psi-sec)"! 
Elastic limit Oyp 150 psi 
Failure by simple: 

o compression 400 psi 

e shear 117 psi 

«tension 200 psi 


These values are not considered to be rate or time 
dependent for this model formulation, although in general 


they are. 
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CHAPTER 3 
Analytical Modelling 


In the existing literature, two basic mathematical 
models were found for the ice sheet - circular pile prob- 
lem. Both models will be quickly reviewed here and the 
features and limitations discussed, 

3.1 Model of Frederking and Gold 

This model, presented by reference [11], draws its 
basis from the earlier mathematical treatise of Noble and 
Hussein [24]. In this earlier work, Noble and Hussein 
derived the exact solution for dual trigonometric (Fourier) 
series that arise in the solution of a mixed boundary value 
problem, such as that one depicted in Figure 3-l. 

The assumptions and limitations inherent in the analysis 
of Frederking and Gold are: 

- Interface friction is not considered 

- Plain strain is assumed 

= An infinite elastic medium is used 
Unfortunately, the series presented as the solution is red- 


ucible only for the special case: 


£, (-›') = (1- 2v) (3.1) 
(orimed quantities refer to properties of the pile) 


Frederking and Gold, however, departed from this analysis early 
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and used the stress function ቀ suitable for the boundary 


conditions; 


ን ф vr TAS ሦነኑ ፍን ፄነላ ርኃ - (1- ፡› e lar + 
+ 21456) D, lae Ae 1 | Dn cosne (3.2) 


n= 2 





[where E is total force on the pile concentrated at 
the center of the rigid pile| 
to derive the radial stress and strain. The infinite series 
of (3.2) was then greatly simplified by assuming a Poisson 
ratio of .5, whereupon the representation for radial strain 


at the interface becomes: 


26€, (a, 0) = po + TE (3.3) 


The constant Dp can be uniquely determined in terms of F. 

Frederking and Gold then go on to develop a refinement 
of (3.3) in which temperature effects, viscoelasticity, and 
a relation between radial strain and penetration rate are 


discussed. Finally an equation relating total force: 
| € E 
&- k(mLc) 8, [E | o^ (3.4) 


where K(T) is a temperature correction function 
v is a constant modifying yield stress by 
strain rate effect 
and L(8)is a geometry function based on (3.3) given 
by: 








Sul 


7 ነተኛ 
| (9) = | (cos Ө) de (3:5) 


here is the "contact" angle 
3.2 Model of Ross, Hanagud, and Sidhu 

Ross, et al., [27] considered adaptation of one of their 
earlier works [34] concerning elastic-plastic plate impact 
study to the problem of an ice floe surrounding a rigid 
inclusion. The model chosen for this study is based on a 
finite difference scheme for the solution of a plane-stress, 
dynamic problem. The problem is formulated with a two 
dimensional cartesian network. 

The method presented contains the following assumptions, 
although it is clear that the model is not restricted to 
these: 

- pile is rigid 

- loading is by edge dislocation or stress on a finite 

plate 

- perfect adhesion of ice to the pile 

- plane stress situation 

The method of solution of this dynamic problem is 
presented in reference [34]. It consists basically of solving 
the continuity equations for the body forces and hence particle 
accelerations. The accelerations are translated into strain 
rates through kinematic relations. Finally the stress-strain 


relations are used to determine stress rates, These stress 
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rates are then integrated numerically into new stress quantities 


whereupon the process is repeated, 


3.3 Criteria for the new model 


For the developement of a model for analysis of the 


circular pile problem, it is felt that the following items 


should 


be incorporated into the ideal: 


The analysis mould present variables to include 
stress, strain, their rates. and displacements 

for the complete spatial distribution and at various 
time intervals. 

The model should at least be a two dimensional 
representation and be capable of plane stress or 
plane strain variation. 

The model should be flexible enough to accept the 
various material representations of ice including 
non-linear and time dependent ones. 

Geometrically, the model should be capable of handling 
mixed boundary values, interface friction, and finite 
or (semi) infinite sheet formulation. 

Fracture or material failure can be easily be incor- 
porated, 

The model must of course be simple to use and econ- 


omical to apply. 


Unfortunately, these ideals are often mutually exclusive, 


In the next section and in the subsequent chapters, a model 


Will be developed that meets at least several of the criteria 
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and will be shown to be useful for application to several 
ice - circular pile interaction problems, 
3.4 Developement of the New Model 

The model to be developed in chapter four is a departure 
from the one described by Ross, et al., in section 3-2. The 
approach used by Ross seemed interesting in that it is capable 
of meeting several of the criteria required. 

The method of solution remained critically dependent on 
the time integration procedure, however, and any numerical 
inaccuracies or instabilities would propagate throughout the 
solution. The dynamic problem also seemed inappropriate for 
study of slowly varying stress fields. Lastly, Ross gives no 
results for his dynamic problem. 

The model proposed in this thesis, although similar in 
conception to Ross', exploits the radial symmetry associated 
with the physical problem to work in only one dimension. 

The method of solution is greatly dissimilar in that it does 
not treat the dynamic problem and uses a linear system of 
equations to arrive at the final result. The new model 

also treats the case of linear viscous behavior and has 


capabilities of extension into other material behavior. 





| 
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CHAPTER 4 


Developement of a Finite Difference 
Elastic Constitutive Model 





h.l General 
The equilibrium equations by Wang, reference [31], 
are applied to the element depicted in Figures 4-1 through 


4-4. In two dimensional polar co-ordinates: 


96, EIS 6) 2ዮ6 br- 6 s ° 
_|_ 6 O Gre £4 E ° 
995 ar + + Ore = Q Ve ae) 


The right hand sides are body forces; the dots indicate 
differentiation with respect to time, 

Since sea ice is a transversely orthotropic material, 
variation of mechanical properties in the Z Co-ordinate 
direction will be due mainly to salinity variation, temper- 
ature profile, and non-uniform thickness. For the purpose 
of this study, however, the ice sheet will be considered 
to be uniform in thickness with zero temperature gradient. 
(Extensions of the present model to account for these effects 
are possible.) For this reason, developement of (4.1) and 
(4.2) into a two dimensional plane stress model is an 
appropriate choice, 

Before proceding, it is perhaps worthwhile to state 


the assumptions of a plane stress formulation: 
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Zelte 


Stresses acting on a small element ABCD in polar 


coordinates 
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FIGURE 4-2 


Variation of direct and shear stress with angular change 


PS ae т 





ZU 





FIGURE 4-3 
Compatibility relationship 
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FIGURE 4-4 


Sign convention 





38 


- CS, Órz,007 are zero'on both faces of the sheet 
- Therefore, Oz, Org, O are zero throughout the sheet, 
- The state of stress can be specified by Dr, De, Dee 
only and variations of these quantities are 
independent of vertical position. 
4,2 Developement of the Plane Stress Elastic Model 
After Wang [31], the strain-displacement relations 


for small deformations are: 


Er = Hur (4,3) 

= Ur | a a4 
ጠን ee Ue (4.4) 
ro + Ur + ŽU- Eë (4.5) 


Stress-strain relationships dictated by the plane stress 


assumptions: 


| ዩኑ = — ( Or - P Os ) (4.6) 
ዩ6 = + ( e-r ô) (4.7) 
Cre = 2(1+p) 6 (4.8) 


E ro 
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Combining (4.1) through (4.8), the following five equations 


are obtained, 


EE s 
Betr Cre =O pco 
>u “ጭ (.6ኑ ”“ ሦዕሬ) =O (4,11) 
т Че + ү ፳6 ሠ “ጭ (69 ”ዕኑ) - ር . (5,32) 
+ 2 ur + 2 Uo- Lu, - BUEH б. О (аз) 


The body forces expressed in (4.1) and (4.2) were 
set to zero since it is not the purpose of this study to 
expose dynamic considerations or impact shock. For the 
range of validity of this model, inertial effects are 
negligible. 

Equations (4.9) through (4.13) (or their cartesian 
counterparts) could be attacked directly in finite difference 
form, but the ensuing number of equations would be quite 
large and difficult to apply. 

If boundary conditions are not to be mixed at a 
continuous boundary, a Fourier simplification is possible; 

Let O,,O5 and Uy be expressed as a Fourier cosine 

(even) series, 


Let Ore, Us be expressed as a Fourier sine (odd) 
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со 

Or(v,9) = > | Or, (Y) cos MO (4,14) 
ነላ = ረጋ 

Og, Uy, are similar 
с© 

Ore (1,8) = > | Ore, (4) 51ከ ኽፎ (4.15) 
ma O 


Us is similar 


It is then possible to express (4.9): 


oo 





zl ፖ, Or, (Fr) cos m6 |+ T | 2 Delt sind | + 
1 253 Os, (P ees me]- 2 [ > 'бе (©) созт o] =o 
(4.16) 
which reduces to: 
x > 206, cos m0 + m» Öre cos mO + 
| + y E (Or, - 66,,) cos mO (4.17) 
Finally: 


)2- Or, + O + -= ( Oy, - Os.) = 0 
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k Or, + — Ore, + + ( Or, - Son cos 8 


| 
3 


E Or, + ዴ Ore, + = ( Or, - T cos 20 
(4.18) 


Thus Or, can be expressed as an equation relating only 
Like wavenumbers. The set of solutions unique to each 
wavenumber can then be reconstituted into a Fourier series 
by using an appropriate number of waves, 

Equations (4.10) through (4.13) can similarly be trans- 


formed into an infinite number of sets of five equations: 


d m m- 4,1 
ፅዮ Orm H Y Oro, т = r Dem mS 
= Oom + = Orem + = Oren 

Ò . 1 " 
Sr Arm E ( Orm ` P Oom) oa О 


an m | 
r Urm tF Uom EE OE ? Or, ) 


m ዕ | 2(1፥ሦ) 
Zy Urm S. Uom TTF Mem 7 — Ore, 
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Ò o m in the fourth row of (4.19) can be explicitly 


expressed as: 


Oe. F- Ur, = Em Uam ү Or (4.20) 


Since no equation contains a derivative of Som this 


variable can be eliminated, yielding: 




















2 6,+( nx JO, + ( + GE ) Ore * ( E-) ur E 2) Ue (4,21) 
غا ءوس‎ ነ6ሬ፥ (E) u 
ፅዮ “ትፄ r RES re y? 9| ዕ 
ur ro rt Hu 
9. -2(1+P) -M -| 
av Let) o JO,+ ( E )Oro +=) Ur A+) Ue 
Let the vector Z represent: O < vx. < cO 
Be 
É = Orom (4.22) 
— JA ü Urm 
don 


| Then (4.21) becomes the set of linear equations: 


»[11Z,+[AlZ,=0 ee 


where [I] is the 4x4 identity matrix and 








— በጨ” am 
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where EN is expressed by the equation: 














EI m E -Em 
r r በ6 r? 
-mp 2 КОЕ -Em 
r Y ES ns 
D si - (1-9) ሦ vin ae 
ج‎ 
-2({+Р) -m -] 
О = ከ፳ 


(4.23) represents an infinite set of 4 simultaneous first 
order linear differential equations. It can be seen that 
two spatial variables would result in an equation relating 
a derivative in © as well. 

The system of equations expressed by (4.23) can only 
be solved numerically. 
4.3 Numerical Solution of the Plane Stress Elastic Model 

The method of finite differences will be exploited to 
solve equation (4,23) 

Considering the finite mesh scheme presented in Figure 
4-5 and in Figure 4-6, the appropriate difference equations 


valid at the nodes are expressed: 


E (32,4 Z-ZÁASNZ.-o 


valid at node 1 (4.25) 





4A 


Urs Ug Or 
Seine 





FIGURE 4-5 


Finite difference representation 


FIGURE 46 


Node numbering scheme 
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mU E Z - aZ + ھا‎ Z-Q (4.26) 


valid at the I-th or last node, 


and the central difference 


AAA 


valid at all interior nodes. 
Recasting equations (4.25), (4.27), and (4.26): 


Equations (4.28), (4.29), and (4.30) are then loaded 
into a suitable simultaneous equation - matrix representation 
as shown in Figure 4-7. 

This representation is then expanded into the full coeff- 
icient matrix shown in Figure 4-8, It should be noted that 


even the simplest scheme, a three node model, will require 
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Expanded matrix representation 
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reduction of a 144 element square matrix (12 equations in 
12 unknowns). Even a numerical scheme such as the Gauss 
elimination or the Gauss-Jordan reduction scheme is too 
laborious for hand calculation, so a machine solution must 
be employed. 

Figure 4-8 represents the expansion of equations (4,28) 
througn (4.30). The elements noted represent the corres- 


ponding position in (4.24). As an example: 


2 


A. r r is evaluated at node 2 


Elements labeled B. or E represent modification to [Al 
required by forward difference and backward difference 
equations (4.25) and (4.26). The matrix is a square matrix 
distinguished by a band twelve elements wide. Even for a 
small number of node equations the ensuing sparse matrix is 
wasteful in space and computation time. Consequently the 
representation in Figure 4-8 is further compressed to a 
band storage array depicted in Figure 4-9, 
4.4 Non-dimensionalization 

In order to produce array elements of approximately 


the same order of magnitude, (4.23) was converted to non- 


dimensional form by the following conversion: 





on stress ° = © (non-dimensional ) 
REF 


where Once = E eer 
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FIGURE 4-9 
Matrix in the band representation 
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E REF is taken as the average Young's modulus through 
the thickness, and length is non-dimensional by 


L = 


=. = u (non-dimensional) 
L ner 


where | = R or inner radius 
Equation (4.23) becomes non-dimensional and the coefficient 


matrix (4.24) becomes: 


























I-P m -| > ሃላ 
RAO RAD RAD? RAD? 
ST 2. - m m 
RAO RAD RAO? RAO? 
D US (4,31) 
D O RAO RAO 
-M = | 
O “2ሯ(1፥ኮ) A EIS 
Here 7 | + E l) — (4.32) 
RAD = Rin h 


and will vary within the array according to node position i. 


4,5 Machine Solution 

Equation (4.23) was implemented in a Fortran language 
Subroutine called elastic. For reasons discussed in the next 
section, the number of nodes used was 400. Consequently the 
core storage requirements in double precision amounts to 


approximately 60K words. 








JL 


Solution of (4.23) was by means of an T.M.S.L. 
subroutine called LEQTIB [16]. This method of solution 
of a system of linear equation is described by [23] as 
the Crout factorization. The Crout factorization is similar 
to Gauss reduction methods but the lower triangular decon- 
E ition is used instead. Partial pivoting and row normal- 
ization are employed. 

LEQTIB was slightly modified into a double precision 
routine and labeled dpband by the author. This routine is 
included as an integral part of subroutine elastic and is 
listed in Appendix A, 

Subroutine elastic was written for the multics system 
implimented on a Honeywell 6180 used by-the Massachusetts 
Institute of Technology Information Processing Center. 
Approximate CPU time for one solution was 16 seconds. 

A flow chart for subroutine elastic is presented in 
Figure 4-10. A symbol table, description of common blocks, 
listing, and a user's guide are included in Appendices A-l 
through A-4, 

4-6 Numerical Stability 

The number of stations chosen is HOO, Table 4-I and 
Figure 4-11 illustrate the convergence of the numerical 
solution as the number of stations is increased for the 
particular case of zero wavenumber. Having checked the 


convergence of the numerical solution, the question now is 





De 


Subroutine Elastic 


Set material constants 


BER 75 A, 5 TEES 
Zero arrays C, D, R | Zerocar || 


Derive fourier coefficients | Loader l | Fourier] 15355 | 
for boundary conditions 


Display fourier coefficients | writzo | 





of boundary conditions 





Read wave numbers desired 


-—- — 
| Load coefficient matrix A y ES 
Load boundary conditions into | Loadbdy | 






coefficient and solution matrices 


Solve | Dpband | 


store wave solutions 


| 
! 
| Calculate stresses, strains, апа | strain || 
1 
f 


present output 


Combine series “ተ፦መ= | 


E E no s wavenumber = wavenumber desired? 


Display total solution | write || 


© FIGURE 4-10 


Block diagram for subroutine Elastic 
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Table 4-I 


Wavenumber Maximum Oscillation 


(per cent) 


О 0.0075 
1 1.7500 
2 0.0000 
3 0.0000 
4 0.0000 
> 0.0000 
6 0.0000 
T 0.0000 
8 0.0000 
9 0.0000 
10 0.0000 
Radius — outer = 20. 
Radius = inner = 1, 
u, = us = O. inner boundary 
Op = -.66666667 outer boundary 
Oro= 0» outer boundary 


Maximum oscillation in numerical solution by wavenumber 


(solution carried to five digits) 
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SE 


whether the numerical scheme is consistent, i.e. whether 
this convergent numerical solution is the true solution. 
| This is done in the following section. 
| 4.7 Comparison with Analytic Solution and other Tests 
| Fortunately, the analytic solutions to the relations 
expressed in equations (4.21) can be obtained and are 
| presented in Appendix B. The analytic solutions and their 
computed counterparts are illustrated in Figures 4-12 
through 4-14 for wave numbers 0, 1, and 2 respectively 
for the case of radial stress at R out/R in = 20, 

Test A This was run to determine the effect of radius 
on the inner stresses, It would be expected that for a 
stress loading of the outer boundary, the importance of 
the radius on the inner boundary would decrease with 
increased plate diameter, This has been determined to be 
so numerically, and some results are presented in Figures 
4-15, 4-16, 4-17, and 4-18. It is also noteworthy that the 
| higher wave numbers have very little effect on the inner 
stress solution. A conclusion here can be drawn that will 
Save much computational effort if only a limited amount of 
data is desired. 

Test B This was run to determine the sensitivity 
of the solution to the Poisson effect. The results for 
four wave numbers are presented as Figure 4-19, Wave 
number one appeared to be almost entirely insensitive, 


whereas wave number zero appeared to be the most sensitive. 
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Numerical and analytical solutions wavenumber zero 
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Numerical and analytical solutions for wavenumber one’ 
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Numerical and analytical solutions for wavenumber two 
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FIGURE 4-19 


Variation of radial stress at the inner boundary 


as a function of poisson ratio for several wavenumbers 
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The next two tests were of load magnitude and shape 
and designed to test for linearity. 

Test C The load magnitude was varied by integral 
| multiples, retaining the same distribution, Results were 
exact multiples as well. No results are presented. 


Test D The test of load shape was more of a test 


| of the Fourier subroutine which converts boundary condi- 





tions into series. The following radial loading was 


supplied: 
A fe #3 
6r 21.cosQ = 1.-1. cos © = 1. 
Ore = O = O. =(), 
Ur inner = O. = С). P 20 
(16 INNER = O, = O. = - 


The results were excellent and are presented in Table 4-II. 
Test E The last test concerned the radius at which 
the finite plate could be considered infinite for practical 
purposes, This test is similar to Test A. This test was 
somewhat inconclusive since above a certain radius some 
solutions became unstable again and innacurate, A practical 
| limit of radial spacing was chosen at n=400, since the 
| marginal accuracy of finer differences was outweighed by 
| increases in computational time, cost, and round off error, 
| Results of this test are presented as Test E in Table 4-III. 
4,8 Comparison with Ross! Solution 


сарин Йан O рыл; 


In Section 3.2 a review of the analysis of Ross was 





Table 4-II 





Radius Run 1 Run 2 Sun Run 3 
1.00 9.72992 -8.2323 1.49762 1.49759 
2.90 5. 34873 —4.2905 1.05823 1.05824 
4.80 3. 36513 -2.3445 1.02063 1.02041 
6.10 2.46312 -1.4538 1.00992 1.00989 
8.60 1.95948 =() 95 361 1.00561 1.00551 
10.50 1.64134 -0.63805 1.00329 1.00328 
12.40 1.42509 -0.42309 1.00200 1.00200 
14.30 1.27047 -0.26917 1.00130 1.00119 
16.20 1.15562 -0.15491 1.00065 1.00065 
18.10 1.06816 -0.067790 1.00037 1.00027 
20.00 1.00002 -0.11921e-04 1.00000 1.00000 
Up, Ug=O dp, ug*O TG 
inner inner inner 
outer outer outer 
OrzcosO Orz1-cosQ Ор= 
Rout=2O 887 2509 2959 20 


Test of élastic routine by variation of load 


shape 
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Table 4-III 
Dr outer -10 
Or inner 5% deviation Instability 
wavenumber @Rout @Rout @Rout 
О >100 > 100 15 
1 2100 50 10 
2 2300 ? 100 [5 
3 25 » 100 7200 
»4 <25 >100 PUT 


Practical extent of finite condition at certain wavenumbers 


by variation of outer radius 
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presented. This paper, reference [27], also included an 
example in which the loading depicted in Figure 4-20 was 
analyzed. Below the stated elastic limit, it was possible 
| to compare the model developed in this chapter. An effort 
was made to correlate this result with a similar loading 
| applied to the method developed in this chapter, The 
comparison is presented as Figures 4-21, 4-22, 4-23, 
| Although the basic shape and critical values of cross- 
| Over are Similar, discrepancies were noted with respect 
| to the hoop stress and also with respect to the order of 
magnitude. By changing the displacement loading magnitude 
"to „Ol inches vice „Ol feet, the numbers were much alike, 


| No explanation was found for the difference in hoop stress, 


but from equation (5.20): 


E E 
De," Tr Arm * ES Uom * P Orm 


Evaluated at Rw : Us, Un = O (boundary conditions) 


| produces 


| Som = P Or, 
| 


kor a non-zero value for any non-zero radial stress, 

x This loading was also run with 3, 5, and 11 wave numbers 
| апа no effect was found on the value of inner face stress, 

| although Significant variations existed toward the edge. 

This is as expected, since the influence of the trigonometric 


components w»| do not substantially contribute to the center 
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U=U ,cose 


Forced Boundary 


FIGURE 4-20 


Discontinuous edge loading of a finite plate about 


a rigid inclusion (after Ross, et al, [27] ) 
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stress as the edge of the plate is extended toward (perceived) 
infinity. 
| Figure 4-24 displays the radial loading represented 
| by three Fourier coefficients and the original loading. It 
| can be seen that even for this discontinuous loading, three 
| wave numbers presents a good approximation. 
| 4,9 Peculiarities and Range of Validity 
This elastic formulation is derived from the basic 
| plane stress equations and consequently should be useful 
where plane stress is appropriato, The model achieves 
Significant accuracy up to a ratio of outer radius to inner 
| radius of 25:1 . Beyond this, the solution diverges from 
the analytical derivation. 
The worst case is at wave number one (Figure 4-25). 
; It is doubtful that the infinite case can be modeled by a 
| finite difference technique in this problem, since a 
x boundary loading for the first wave number at all radii 
Will produce a significant solution at the pile, Consequently, 
this elastic model is useful for radial limits of outer radius/ 
| inner radius ratios less than 25. 
| The buckling limit of a thin plate is an important 
| consideration. The analytic derivation of one buckling load 
is presented as Appendix C. The equations of plane stress 
do not allow for deflections in the third dimension, therefore 


| they will not predict the failure by this mode, External means 
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FIGURE 4-24 


Displacement loading (radial) applied to boundary 


of Ross! solution - comparison with first three 


wavenumbers 
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must be applied to the program to preclude this occurance. 






The plate is assumed to be attached to the rigid in- 


! 
1 


clusion (frozen-in problem). The partial movement of either 


| 
D 


| 
| 
| 
| boundary cannot be supported by this method, although it is 


| conceivable that a more general two-dimensional model could, 





Finally, the model can support an application of forces 
or displacement at either boundary, which can in turn be 
discontinuous. The degree of refinement is limited by the 
| wavenumbers considered, The accuracy is limited by the 


| circumferential spacing chosen, 
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CHAPTER 5 
Viscous Extension of the Model 





| 

1 
| 5,1 General 
Considering the simple rheological model in Figure 5-1, 
| the following differential relations hold: 


ር = С, viscous) + C, ELASTIC (5.1) 
Ce = ር. viscous + a ELASTIC (5.2) 
Gs = ዩኤ viscous + Er ELASTIC (553) 

where C viscoos = rm ( 6) (5.4) 


The Maxwell model is chosen for its simplicity and 
` because it is possible to generalize into a non-linear 
creep law. Additionally, some investigators have modeled 
| sea ice behavior this way experimentally, and physical 
values are available [33]. 
| 5.2 Incorporation of the Viscosity Relations 
The stress-strain relations of (4.6), (4.7), and (4.8) 
can be expressed using (5.1), (5.2), and (5.3): 


É, = = (Ó, - Gei + = (Š, - ሦስ) (5.5) 


° e 


(= <(09- POr) + (O, - Or) (5.6) 
( _ 2% (| + >) Öre (5.7) 


e E 








E 





viscouy 
HE 


LARI Na 


Ly relaxation of 
stress 


/ / p.d.-permanent 


E < 


deformation 
/ T.s.-recovered 
RE — X) / strain 
| Zei 
ሽን s -5-' 


GU RE 5-1 
Model for Maxwell liquid (after |8] ) 
Al=Als Al, 
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These relations can be incorporated into a series of equa- 


tions as in section 4.2 yielding: 


a§2(t]Z, + (AlZ,3 =[EIZ. os 


| Where: O < w < co 





— 


Al is the same matrix as expressed by (4.24) 








` And 
k E. ch 
-t -E ma 
O © 
r? p? 








[E] _ Š S S S (5.9) 


ፎል (1+ሠነ- ቼጄ ክለ (1ት») 
v ^". 


О O 


O det (1+p) o 


O 
E 


5.3 Solution 

Numerical solution of the series of differential 
equations (5.8) is carried as before, only the unknown 
variable, SÉ » represents rate at time { . The value 
| ZÉ used e right hand side of (5.8) is the elastic 
| Solution for É . Note that boundary conditions of load 


— M 
rate must be employed. The solution yields: 


ol A CONSTANT 
at | Z í = Cw 9) 








9 


| This differential equation in time can be solved using the 
Euler or "one step" method to arrive at a solution for the 


next time step, where the process is repeated, 


2 Т t+At 


(Zat + Z. = گر‎ (5.12) 


4 


Non-dimensionalization with respect to another variable 
is required. The variable chosed was ol, the viscosity 
coefficient, with units of | T М. 

If Xreference is chosen, and mass is preset by consideration 
of stress non-dimensionalization (Section 4.4), then derived 


time: 


| ° 
We к ы е (5.12) 
E E ref OC REF 


Therefore any value of nonedimensional time is related: 


T T ~T puysica 
Teer 


As an example, using typical values for sea ice [ ], 


(5:13) 


pee = 2083 psf 
Xo = | 27 -3 B 
eec “= 1l.2788x10O (pst - sec) 
implies E: 
E = | 
and x= |. Thus “lage = 2007 sec 


or one unit of nonedimensional time 
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A flow chart of the numerical computer solution and 
time integration is presented as Figure 5-2, 
5.4 Convergence in Time 
| The Euler method chosen for time integration assures 
| convergence if the time step is small enough, Since it 
| is not possible to analytically determine the solution, as 
| was done for the elastic case, one approach is to run a 
| general solution for increasingly small time steps until 
[ no numerical improvement is achieved, Unfortunately the 
solution time can be prohibitive. The results of two time 
convergence test series, for wave number zero and one, were 
found to converge at a time step of .001 T » beyond which 
" significant improvement was not achieved, The time increment 
to achieve this numerical stability comes out to be a constant 
x independent of stress, for the upper bound prediction, 
5.5 Range of Validity 

The viscoelastic formulation presented in this chapter 
|. is not valid where the elastic theory is inappropriate or 
where the numerical limitations of subroutine elastic do 
| not apply. The range of validity is further restricted 
x on how well and on which range of parameters ice can be 


modeled as a linear viscoelastic material, 
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Subroutine Visco 





Read time step 


Read stop time 


Read time solutions to 


be displayed 


Save elastic solution 

Load modified solution |Loade | 
vector 

Load rate boundary conditions 4 


solve system of equations 222227 
for rates 


solve differential equation 





Increment time step 






> display interval reached? 
Display solution [===] 
> time=stop time? 


SG ON ን ጋ 


Block diagram for subroutine Visco 
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CHAPTER 6 


Analysis of Ice - Circular Pile Behavior with the Model 


The previous chapters have dealt with the need for an 
analytical tool and the construction of a mathematical 
model, Without exhaustively analyzing the many combinations 
of loadings and important parameters that affect the circular 
pile problem, some interesting results will now be developed. 
6.1 Material Behavior 

Katona and Vaudrey [18] have presented a logical seq- 
uence of obtaining ice material behavior experimentally for 
a particular geometric configuration. The results of these 
data are then fitted into curves which can then yield a range 
of validity of ice material law as discussed in Section 2.3 
and with Figure 2-2, 

Of course, since this model formulated in this thesis 
responds only viscoelastically, the transition from visco- 
elastic to viscoplastic behavior cannot be obtained. The 
model can, however, define the elastic-viscoelastic transition, 

As an example, consider the loading situation depicted 
in Figure 6-1. Physically, this could represent the rigid, 
frozen-in pile loaded by a large, thermally expanding ice 
sheet, Neglecting any temperature effects to viscosity and 
elastic modulus, the loading curves of Figure 6-2 are con- 
structed. Stress-strain diagrams similar to Figure 6-3 show 


the viscoelastic creep behavior. The apparent relaxion of 
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It 
SESS 


Viscosity coefficient O< 


11 
= x 


Elastic Modulus E 
Non-dimensional quantities, 
در‎ 


see text in Sections 4-4 and 





FIGURE 6-1 


Circular ice sheet — rigid pile geometry and loading 


configuration used in viscoelastic behavior analysis 


of Section 6.1. 
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Effective 8 
Strain ОШО) 


E 2 3 
Cen- (€ s 6.64 q^ 


100 Elastic modulus = |. (nondimensional) 


3 


Viscosity coefficient - E (nondimensional ) 
Variation of strain at pile-ice 
interface 


50 Loading condition of Figure 6-1 applies 


Outer boundary conditions 


| 20 
Radial load applied 
| @ = Bor 
10 
ቱጌ. 
2 
= _, 2232107? 
2 
=_.111х107? 
1 
О sch ac 3 en $5 
Load duration - non-dimensional time 


FIGURE 6-2 


Viscoelastic behavior - Effective strain VS. load duration 
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Effective 3 
Stress See ki 
MRAZ 2 2 
Bett "Lal De" Dr Da 365) 
2 
Non-dimensional time 
OFA Ne „О 
$15 
cT 
„2 
О 
О eL i2 > ፅፋ SCH 
Ee 1. 
2 „2 2 2 
FIGURE 6-3 Ceff =(Er+* ES* Ep Est Eno) 


Viscoelastic behavior — stress-strain diagram showing 
creep for loading condition of Figure 6-1, at the pile- 


ice sheet interface. 
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elastic modulus is then limited to a value of 10% decrease, 
This then becomes the transition boundary between visco- 
elastic and elastic behavior of Figure 2-2, 

In the example, loadings of those shown in Figure 6-2 
are used to define a region of validity in Figure 6-4, 

Using actual values such as those found in Section 
2-4, it can be seen that a time of approximately 10 minutes 
defines the region at low stress levels (below yield). 

The effective stress and strain plotted in Figure 6-3 
are defined: 2 | 2 2 2 
O EFF == (0, + Og - OF. G4 +3 Ore) (6.1) 


2 2 2 
Ess Lë + be + & 6, + bre) (6.2) 


For this loading condition, it can be seen that a 
value for modulus of elasticity is about one half the uni- 
axial value. This will vary as well along each point in 
the radius of the ice sheet, Interestingly enough, the 
relaxation to 90% of its former value occurred at the same 
time, despite the variation. Under a more general loading, 
not confined to the zero wavenumber, it is likely that the 
material will behave not in the same region of Figure 2-2 
Over its entire geometry. 

6.2 Under Rate Loading 

The example in the previous Section was subjected to 

rate conditions on the loading. This was done to test the 


capacity for this parameter in the model, The load-time 
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Load m 
Magnitude х10 


20 | 
| 

| 

| 

LO ° 
2 ቁ 

Elastic Region Viscoelastic Region 
2 Y 
2 
SL 
О » d. „2 each . 1 5 


Load duration — non-dimensional time 


FIGURE 6-4 


Viscoelastic behavior — Elastic-Viscoelastic transition 
from load curves of Figure 6-2. By definition, relaxation 


of elastic modulus is limited to 10%. 











88 


curves for the loading of Figure 6-1 are presented in Figure 
6-5 for three loading rates. The greatest rate is very 
close to the maximum of 5000 psi/min suggested by Katona 
and Vaudrey to avoid dynamic effects, No general conclusions 
are drawn here since the size of the initial load and the 
various rates have little significance without a specific 
physical problem to simulate. 
6.3 Analysis of a Friction Coefficient 

The loading situation in Figure 6-6 was applied to the 
model to determine the ratio of shear stress to (compressive) 
radial stress at the rigid pile interface. This could be 
construed to be an indication of a static type of frictional 
coefficient, but since the model also assumes capability 
of the interface to support a tension, in this region "fric- 
tion coefficient" would be meaningless, 

The analysis proceeded with those values that this edge 
loading produced, Defining the coefficient of static friction 


as; 
Ore 


Ms = Mee (6.3) 





it can be seen that this coefficient will vary about the 
interface at least to the point at o: T where Oy becomes 
tensile and a coefficient of friction becomes meaningless 
(Figure 6-7). 


At a value of Oz .31 radian (given this loading), 


Ms = .5 
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Effective 3 
Strain x10 1 
22 22 
Ceff = (Er*€e*€r fo" Ero? 
re Load a 
-20.x107^/T 
24 
E -2. 307” /ጥ 
: 0. /T 
O T 2 vo ES 
Load duration — time 
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FIGURE 6-95 
Visco-elastic behavior — Ceff VS» time curve under 
three loading rates, at the pile - ice sheet interface. 


Loading condition of Figure 6-1 applies. with C -.1111x1077 
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Rout = 20 * 
Op =--333 nase outer boundary 


Ong=-- 333 sin® cos@ outer boundary 
Tor TESST 


and O. otherwise 





at the inner boundary or pile-ice 


sheet interface: 


Rin = dl og 

Ur =U, ze Om considering "friction" 
and 

Org =Up -O. "frictionless" 


FIGURE 6-G 


Loading condition circular pile - ice sheet interaction 


for cases discussed in Sections 6-3 and 6-4. 
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Ratio 


1.5 / \ 


1.0 


e A 


p 


| compression tension 


Ome 020,912, 1,25773.571.1.885. 2.199 2.513 2,827 3.1 


Angular position in radians 


a 


frictional considerations valid 


in this region 


FIGURE O7 


Ratio of shear stress to radial stress at the circular 


pile =- ice sheet interface, loading of Figure 6-6 applies. 
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Now if there is no ice - pile adhesion present 
and if the frictional coefficient of smooth sea ice on 
a smooth steel circular pile is about „3 as some inves- 
tigators indicate, then it is clear that for regions 
beyond .3l1 radian Ug will be nonzero. A situation in- 
volving mixed boundaries and separation will than be evident. 


By changing the loading to: 


6, = - O. COS Ө 
Ore z- Ore Sin Ө 


it was found that the coefficient of friction at the inter- 
face changed only slightly. Since for the wavenumber zero, 
only normal stresses are conveyed (being an odd function) 
and the higher wavenumbers lose their influence with increased 
radius, frictional forces are important only in the region 
| le | < .31 radian for non-axisymmetric loads. 

A conclusion can be drawn that for the model of the 
Circular pile - ice sheet interaction and unless ice adhesion 
is present, frictional behavior is confined to a small area 


preceding the pile and may be unimportant in analysis. 


6.4 Analysis of the Behavior of Ice Floe =- Pile Interaction 
Comparing Adhesive and Non-adhesive Behavior 


Since the influence of friction at the circular pile is 
Minimal, it is necessary now to compare the behavior at the 
interface using complete adhesion (frozen-in situation) and 


and using complete circumferential freedom, 








SE 


The loading of Figure 6-6 again was applied to the 
model to consider the effect on the radial stresses of a 
"frictionless" boundary. Since the loading is the same 
for both, the net force felt by the pile is the same in 
both cases. Figure 6-8 presents the variation of radial 
stresses in both cases, The point where adhesion is 
required is precessed to 1.48 radians. 

Noble and Hussein [24] obtain a value of 1.24 radians 
analytically in their formulation. Severe differences do 
exist in the problem as their model does not support tension. 
Therefore it would be expected that in the new model, 
the radial stress on the compression side would be supporting 
the lost force on the tension side, with basically the same 
shape in the loading curve. 

At any rate, a conclusion can be drawn that the friction- 
less case is the worst as far as local normal pressures on 
the pile structure are concerned, and that the maximum force 
can be as high as three or four times the frozen-in case, 

Of course, the total lateral load remains the same. 

6.5 Locations of Maximum Stress and Strain 

As a final use of the model, information from the loading of 
Figure 6-6 will be presented. Physically, this edge loading 
might be produced by a wave driving force impinging on the 
circular floe, No failure mechanism is included, since 


Some investigators feel that failure is ductile, caused by a 
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FIGURE 6-8 
Comparison of normal (radial) stress &t pile - ice sheet 
interface, with and without friction. Loading of Figure 


6—6 appl ie Se 





Es 


limiting strain, and no reliable theory is available. 
Therefore no estimation of a maximum total force will be 
offered, but a good representation of stress and strain 
distribution can be presented, 

Figure 6-9 depicts the variation of normal (radial) 
and shear stresses throughout the ice sheet and at selected 
angles © . The greatest radial stress is at Q=0° as 
would be expected. The crossover to tension occurs at 6=90°. 
The maximum shear stress occurs at about 54° again at the 
inner face, 

Figures 6-10 and 6-11 show the variation in vertical and 
hoop strain around the pile. If tension strain failure is 
the failure mechanism, it can be seen that failure will occur 
in the sheet vertically (a horizontal crack) and horizontally 
(vertical crack) at angles of 0° and 54° respectively. Shear 
strain follows the shear stress and in this case is 2/3 the 
Org value. Shear strain is shown in Figure 6-12. The values 
of these strains are plotted in Figures 6-10, 6-11, 6-12 only 
at the interface, since this is the location of their maxima. 
It should be pointed out that the basic shape is maintained 
some way into the ice sheet, so that a crack would be oos 
lated at the ice - pile interface and propagate radially 
outward, 

It would seem heuristically that the cracks at 54° 
would occur first due to the greater magnitude of the strain. 


It is then possible that the crack at 0° will then follow due 


the redistribution of stress. Both types have been reported [15]. 
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FIGURE 6-9 
Elastic behavior — Loading of Figure 6-6 — Variation 


of normal (radial) and shear stresses throughout the 


ice plate 
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Elastic behavior - maximum vertical strain location on ice = pile 
interface. Plot of vertical strain vs. angular position at radius 


equal to 1. Loading of Figure 6-6 applies. 
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Elastic behavior — maximum hoop strain location on ice - piie 
interface. Plot of hoop strain vs. angular position at radius 
equal to 1. Loading of Figure 6-6 applies. 
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| Elastic behavior — maximum shear strain location on ice - pile 
interface. Plot of shear strain vs. angular position at radius 


equal to 1. Loading of Figure 6-6 applies. 
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CHAPTER 7 


Conclusions and Recommendations 


El Conclusions 

The model developed in the preceeding chapters is 
applicable to a circular elastic or viscoelastic sheet 
surrounding a rigid circular inclusion. The conditions of 
plane stress are assumed, although plane strain could be 
very easily implemented, Loading can be by displacement 
or stress at either boundary, although mixed conditions 
at a particular boundary are not allowed. The model is 
favorably disposed for inclusion of failure criteria, but 
none are included or proposed. Tension is required at the 
interface to support a valid operation and this corresponds 
to a frozen-in situation. This model can accept time dep- 
endent loadings below dynamic significance. 

The model is very thrifty with computational resources 
and produces a very high degree of correlation with existing 
analytical solutions. The model is capable of considerable 
expansion. 

Some results with a fictional loading predict frictional 
effects are present only at a small region immediately pre- 
ceding the pile. If angular movement is allowed, the resul- 
tant loss of shear requires an increase in normal stresses, 
If a separation of boundary is to occur, and occurs where 


tension cannot be supported, then the contact angle will be 


90°. 
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The range of time where ice can be considered elastic 


(below yield levels) with a uniaxial load is: 


Other loadings could be analyzed. Rate of loading is 
important insofar as time to failure and required force, 
Finally, viscoelastic behavior predicts ice failure (and 
maximum forces) for loading well below immediate fracture. 

Lastly, it is possible to qualitatively see failure 
by limiting strain or shear at locations in the ice sheet 
that are physically relevant. Without knowing the values 
of biaxial strengths of the ice material, it is not possible 
to predict a maximum force before fracture, but it is likely 
that a limiting strain or shear strength precedes ultimate 
crush failure, 
1.2 Recommendations for Future Work 

The model developed in this thesis is extremely rudi- 
mentary and can be used for pile - ice sheet interaction 
under only the most generous conditions, Improvements can 
be made, however, to increase the usefulness into problems 
with more physical application. These improvements by order 
of importance are; 

(1) Driving Forces 
Ice sheets and ice floes are driven by the effects of wind 
and current shears, It would be desireable to develope an 


analytical or numerical procedure whereby these shears can 





| 


— 
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be approximated into edge loading conditions. The size of 
the ice sheet considered in this model would require negli- 
gible wind or current shear itself, so that the two-dimensional 
approximation could remain valid. 

(2) Failure 
The failure mechanism or mechanisms can be easily applied to 
the model if a good choice could be found. To date, inves- 
tigators have proposed Tresca or Von Mises type criteria, 
failure by limiting strain, rate dependent crush strength, 
among others. The behavior of ice is poorly understood in 
biaxial loading and experimental work needs to be carried 
out to properly extend the model into realistic and workable 
significance, 

(3) Two-dimensional Extension 
The model should be given another co-ordinate in the plane 
such that more general geometries can be analyzed. Of course 
this could mean larger quantities of data to be handled. and 
an attendant loss of simplicity. The problem of the mixed 
boundary conditions and different shapes of the ice floes or 
structures could be adjusted to fit physical requirements. 
with a high degree of confidence, 

(4) Multiple Structure 
The model, once in a more general two-dimensional state, 
could be extended to study the effect of multiple structures 
such as support pilings on a pier or multilegged towers, 


Possible reductions in forces by interference could be shown. 
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(5) Other Material Behavior 
Once in a two-dimensional form, the model could be extended 
into non-linear and non-elastic behavior in a straightforward 
manner, The one-dimensional form presented in this thesis 
is valid only where superposition is valid. 

(6) Other Structural Behavior 
Experimental work with the model could be extended into 
non-rigid analysis for study of structural responses and 
forced vibrations. 

Data on the material properties of ice in a biaxial 
state of stress need to be acquired for realistic imple- 
mentation in the model, Once this and the six stated improve- 
ments of the model are implemented, the behavior of the ice 
sheet-structure problem can be fully analyzed and more general 
design equations or curves offered, Additionally, the 
influence of other important parameters such as temperature 
profile or variation of ice thickness can then be included 


or discarded, 





| 


5. 
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APPENDIX A-1 


Program Variables 





Mariables in Common Purpose 
| n number of radial spaces 
| пп number of circumferential spaces 
nnpls nn+l 
ind indicator of boundary conditions 
used (see listing under load- 
bay) 
m wave number 
ru poisson ratio 
el Young's modulus 
pi 3.1415... 
alpha viscous coefficient 
mdes number of wave numbers to be 
processed 
rout outer radius 
rin inner radius 
| idsp 4xN+number of radial display 
| intervals 
ibar termination variable 





| Variables Consistently Used 


irowmax 


number of simultanecus equations 


to be solved 





tstop 





tstep 
time 


disint 


streslim 
strnlim 
rad 


hoop 


A 
Zaff m am 


108 


stop time limit (visco elastic 
solution) 

time step interval 

time since start (elastic: time = 0) 
time display interval of viscoelastic 
solutions 

elastic limit of stress 

limiting strain 

radius 


hoop stress 
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APPENDIX A-2 


Description of Arrays 





‚ Array subroutine 


. /blokl/a(1601,17) elastic, zerabx, loada, strain, 
double precision loadbdy, visco, dpband 


purpose: holds coefficient matrix in band compressed 
form - after call to dpband holds the factored matrix 
/b10k2/b (1601 51) elastic, zerabx, loadbdy, 
visco, loade, boundb, strain, 
dpband 
purpose: solution vector of the system of linear 
equations represented by a 
/blok3/e (204,21) elastic, zercdr, recomb, 
writr2 
purpose: stores total recombined solutions at every 
eighth node (row) and all circumferential positions 
| /oloklt/à (11,8) elastic, zercdr, fourier, 
| loadbdy, writrO 
purpose: stores fourier coefficients of boundary 
conditions for O through 10 wave numbers (row), inner 
conditions (first four columns), outer conditions 
(second four columns) as determined by variable ind 
r1(8,21) loader, fourier 
purpose: stores boundary information from main before 


fourier decomposition 








I 
ኣ 
} 


| 
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ws (21) fourier, simps 
purpose: working vector to be integrated by simps 

/o10k6/x1 (14436) elastic, zerabx, dpband, visco 

double precision 
purpose: working storage used by dpband 

/o10k8/r (204,11) elastic, zercdr, recomb 
purpose: stores elastic wave number solutions for every 
eighth node 

/blok9/t (51,11) elastic,zercdr,recomb 
purpose: stores hoop stress wave number solutions for 
every eighth node (similar to /blok8/) 

/blok10/tt (51,21) elastic, zercdr, recomb, writr2 
purpose: stores total recombined elastic hoop stress 
solution (similar to /blok3/) 

/olok11/rr(1604,1) visco, loade 

double precision 


purpose: working vector 
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APPENDIX A-4 


User Notes 


Program Operation 
The program flow is described in Figures 4-10 and 5-2, 


The material constants E,w, Al are set in either the main 
calling program or subroutine elastic prior to compilation, 
In the listing presented in A-3, E and A are non-dimensional 
with a value of 1,0. The program is not restricted to 
nonedimensional solution, however, and any consistent values 
could be used. 

User inputs are via keyboard and consist of 3 inputs. 

For elastic: 

1) mdes - after reviewing fourier coefficients of 
boundary conditions, the user is asked for the 
number of fourier coefficients to be used in the 
solution. This feature is incorporated to save 
calculation of meaningless higher wave number 
solutions, 

2) The user is queried on whether a viscoelastic 
solution is required. This feature is used to 
permit an elastic solution only. 

3) Display of theta information - This feature 
automatically combines elastic wave solutions and 
calculates circumferential variation for intervals 


of .05 JT] and radial variation of 50 intervals. 
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Display here is set for .OITT and radial variation 
of ten intervals. 


| If a viscoelastic solution is desired, additional 
| 


i 


| keyboard inputs ares: 


1) tstop = desired stop time, entered in units of 
non-dimensional time, 

2) tstep - interval calculation time - This is 
important as a small step must be chosen to ensure 
convergence, 

3) disint - time display interval of wave number 
solutions - It is anticipated that not all small 
discrete time intervals will require display. 
disint should be an integral multiple of tstep, 

Loading Input - rout must be set in the main program, The 
vectors putl and put4 contain boundary condition loading 
information as determined by ind, also sent from main, 
putl - put4 contain diecrete data for loading 
circumferentially through 9=T1T „ (The other half 

is identical by symmetry). 

| | Inner Boundary Conditions Outer Boundary Conditions 


| 394. putl — put? put риш 
| 1 Ur uo Uy uo 

2 Gr Ore б, Cre 

3 ES Ug O, Ore 

4 O, Ore Ur Ue 


Caution must be exercised with ind = 2, since a rigid body 
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translation may be produced. 


Output of Program - Three distinct listings are produced: 


1) 


2) 


3) 





Fourier coefficients of boundary conditions - all 
wave numbers 
Wave number solutions of 6, Oso , Me, Le om Ocg; 
Cere » node, radius for upto 400 nodes. In the 
listing of Appendix A-3, only ten radial intervals 
are displayed. idsp will vary display points., 
Wave number solutions are displayed for the elastic 
case and at each time display intervals determined 
by disint. The headings "stress" and "strain" 
contain an effective stress and strain as defined 
RU I 2 a2 2.5 
6 “፳፳ (6. + бе - бб + 3 Oro) 


E= (E + fa * C fot éro)? 


Stresses by angular position - contains up to 
fifty radial intervals with the total recombined 
fourier series representetion of the elastic 


solution displayed at intervals of .OITT. 


| Other Notes 

Viscoelastic boundary conditions (rates) are preset 
to zero in this program, but are not necessarily restricted 
to this value or even a constant value, Rate information 
must be in size with integration of time step and convergence 


must be ensured overall. 





| 
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streslim and strainlim are preset constants corresponding 
to an elastic limit and a strain failure limit. The solution 
is deemed inappropriate beyond these levels and a solution 
is terminated for that particular wave number. 

The accuracy of the fourier subroutine is dependent on 
the circumferential interval. For complex loading it is 
conceivable that this interval can be decreased by increasing 
nn and of course pertinent array sizes. 

Storage of solutions in the B, C, and R arrays for a 
particular variable are in every fourth row, e,g.; 


solution 


Z =blc,1), bli+1,1), bli+2,1), b(i+3,1) 
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APPENDIX B 


Analytic Elastic Solution for the Case of a Circular Ring 

Wavenumber m=0 

The stress function is known to satisfy the biharmonic relation 
y> © | (B.1) 

In à two-dimensional elastic material. For the case of no 


circumferential variation (as in wavenumber 0), (B.1) reduces 





to 2. 

(2,123 2X 4223) «o 

r Y or 

ae 209 C o ቋቅ 8,2 
Or de* "v de? mar Par (3.2) 
Let f= log r 

dt deat vat 
Then dr dide var 

de, (sty 4 (att at 

v deae] iat? di 

att ot eat dit جل ور‎ 

Sgt m)‏ )د 

LOOSE Fa 

dr as oet IPTE -5ሻ) 


(B.2) becomes 


a^ 4ፋ ME 
iP OM 4ነ" 


The general solution of which is 


ጊነ 
tote ine + Í Cs + CA 


ale gy 


a Belle Te Se 





5 


Finally b= Cir" loge + Cz loar + C; r, C3 loar + Ca 
Now SE = = En (1+ 2 1oar) tC, E = (B.3) 
2 
Oo = Fa tc, )3+ وا2‎ (+2 - Y (8.1) 
Oro * O 
Strain relations 
du LL ሃ 
beige AF wit 
In plane stress 
6. = O 
du | 
RH" pl6, - Gel (5.2) 
~?¢(0,-»6,) (B.6) 
Integrating (B.5) 
Eu = ነ(5.3)ው -* !(5.5)#» 
Into Eu * C, | e (1-3v) + 2. (1-»)(r logr-v)] + 2ር3 (1-ሦ)ኑ 
-C, +») + Ce (5207) 


|(B.6) becomes 


Eu > ር ኒኑ (3-») +2 rv (1-1) loge | + 2C, (1-»)r 


ai tuu s 


(B.8) 


"For (B.7) to have the same representation as (B.8) 


x C. Cs > O 


Therefore, for wavenumber O 


Eu = 2. Cà (\-»)r- C, (48) = 


E 
Orc, > - 


бы 2 Er =- 


d 





(B.9) 


(B.10) 





Wavenumber m=1 154 

From the general solution for the stress function in circular 
coordinates presented by Timoshenko [29], that portion Sa 
responding to the first wavenumber is: 


a - 
T: » rO su. Ó + (kr3 + alr" + o; vloqy )cos € 


E 


| 23 ! ot 
Since O=* dx *¥? 362 


-6 


2 
Go = 2 


Qu) 


Ki 
_ д /! d4 
Gro * 2 X 


Then 6. * (la ii + 2b, r - 2a, | cos O 


Os = (6b r NO r b eos o 


l -3 t ei : 
6,a < (2b, = 2a, + +6, v le, e 
Employing stress-strain (4.11) and integrating with respect to r 


Beet ci 
+ b, log 3 ዩ cos A $ (e) (B.11) 


Employing (4.12) and integrating with respect to © 


Euo x L- (o. * o.) lea" + Sb +O, ር b ከከ 7 
ate ኳ logy + (a, +b, )- b, r =a, >] su Q 
- (ና(5)49 + FC) (B.12) 


< 


From (B.11) = A ` - | (a+ br +0, lag к 
ዓናን r - ars b. loge | sw G + LM 
(8.13) 
From (B.12) E Sue AL E + 10bir -2aır | 596 + 


2e [yr 28e-2a e? ] sie +F' 
(B.14) 
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And also from (B.12) 


-55.. -|- (a,4 ët + Sbr +a, و شه طب‎ 
+P[-b, er (a4 Ee Sor - ዉነ የ”! su. © 
+ уге EC (В.15) 

v 


Y 
Equating (4.8) and (4.5) and using (B.13), (B.14), and (B.15) 


$Q p, p - [46 + (1-4, ]su0/r 


or (ren, Lëteide + (ve'-F) = [Ab (1-H) a] sue 
(5.16) 


Bor) must necessarily be zero 


LE. Eso => Felr 


Differentiating with respect to Ө 
í", í - [46,3 (1-»)a, [cos O (B.17) 
The solution to the differential equation (B.17) 


Gre Assa + Beos o + 120/+ Ea, Jo su e 





Wavenumber > 2 
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From [29] 
“+2 E " 
+ = (Anr"+b, r E + Ee E awe 
ኤሩ 2 zu L 
б, E -Lnlu-1)a,, lr a + w(wal)a, v. + 


+ (Wan-2) bu ro | ces 6 


De Š In(u-N)aur + (n wel), rs n(n Nay, GE 
+ (w-1)(u-2) De v^ сс л Ө 


wm Ms d 


6,e = Lnkn-Uäat? a nus) bur = n (nata r 
-n(u-t) D. የ | Stu ud 


Again employing (4.11) and integrating 


E Ur = - naur” + (u-2) b. клы Ak z: (ur?) b. pu 
‚созуб -е[ калу“ + (ur?) bu ees A UMS 
- (w-2) bu "| es we + £co) 


(8,385) 
Employing (4.12) and integrating 


E መላ” “us 
E ug = \naur E (nsa) bn Na à x (w-A)e v 1 
፡ ዊጨ ኢመ ራሥ ቢላ ሷዚኛ' s u by + an Y 
rnb, et] sane - C&to)ae « FC) 


atu ol 


x (B.19) 
| E ou "à „u | 
x From (B.18) E — = I: [nann tulun- bn rwa ulna bn: 


E > W 2,1 eW 
"ትዮ | о |н о, ү“ ki nlu?) kur- u An Y 


- w(w-3) e, v^ 12 SU. uO + $' 
Le 
(8,20) 


a-a 1 - 
From (B.19) z [ n(u-tlayr + (wet waa) ዉየ - nlus) An RUM 2 
- (u-u-4) b. e" ፋኑ In (n-Naur" v (wa) bu: 
“a nulua - አርን ሀራ r Jese ы 


(B.21) 





10% 


Also from (B.19) 
O W-2 Y 1 “ካ=2 ‘a 
p2., - | nour СОУЛ ` $ (w-4)6, v"]- 


4-2 ሐ I —и-1 
SU nO DAT +ub, v EN A y Y 


iw oO suus OFT -E (B.22) 
| Equating (4.8) and (4.9) and using (B.20), (B.21), & (B.22) 


| F 
t Ns = 
= +E 4 ኑ (dao = Ze o 


ç! + (feodo = E ep e CONSTAN T 


F = C+ Dr 
{- ልፄ5ጨሬ6ት Beos © 


These general solutions must be employed with the correct 
boundary condittons to find the unique constants associated 


with each wavenumber, 
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APPENDIX C 


Buckling of a Circular Ring Subject to Stresses Uniformly 


Distributed Around the Edge 


The loading applied to the condition depicted in Figure 


Cel is analyzed for the first buckling mode. Assuming the 


deflection surface is a surface of revolution, Timoshenko 


[28] gives the required differential equation: 


deb dd Qr? 
CUN MIA ሽማ... 


where Q is the shearing force/length 
and D is the flexural rigidity 


ርን ። Ne sin $ = N. $ 


Ne ^ 


Let = ol 


(C.3) becomes 





2 
с dE + (ar) = O 


Ar 
Lf u= xr (variable change) 
2 
2:15 el $ 2 + 
Then ዚዜ SS, + wT ا‎ 


The general solution of which is 


Applying boundary conditions 
db (rsa) = О clamped inner edge 


dé p 5) ጋ r=b free outer edge 
dv v 


(c.1) 


(C.2) 


(C.3) 


(0.5) 


(0.5) 


(c.6) 


AO 
(055) 
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Figure Cel 


Buckling of a thin circular plate around 2 rigid 


inclusion 
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Since as = Aix ፕ (.) + A, Yu) (0.9) 
The boundary conditions (C.7), (C.8) become 
A, J, (xa) + A, Y (ea) = © (с.10) 


(A, T, (b) + Az Y/ (ab) + 
+ E(A,3,(xb) + AY («b))= 0 (533) 





T, (xa) 
.10) is z- : 05 
i Ba D እየ (xa) | | 
J, (xa 
And (C.11) < (3, (et) - li L. a Yt ul +E( (J, (xb) - (0.13) 
S J (xa) ኔነነ ። 
Y, (o) qu 
Ewe, т (0.14) 


(c.12), (c.13), (6.14) are combined into a single equation. 


The first zero appears at about a value of .42 


a l H3 
5 ' “== 





To apply ጪ። 1, 20 ነ 
Where «47.42/0. then Neg 7.160 
3 
= EN | the flexural rigidity 
12 (1-9) | 


Yields Geg = = z oes EW /ዶ 
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